--- title: "Mt. Vesuvius Lives" subtitle: "INFO 526 - Summer 2025 - Final Project" author: - name: "Seismic Minds - Cedric Destin" affiliations: - name: "School of Information, University of Arizona" description: "Exploratory Analysis of Seismic Activity: Mt. Vesuvius" format: html: code-tools: true code-overflow: wrap code-line-numbers: true embed-resources: true editor: visual code-annotations: hover execute: warning: false echo: false render-on-save: true --- ```{r include=FALSE} # Load necessary libraries options (warn= - 1 )if (! require ("pacman" )) install.packages ("pacman" )# use this line for installing/loading :: p_load (tidyverse,``` ```{r} :: theme_set (ggplot2:: theme_light (base_size = 12 ))options (width = 65 ):: opts_chunk$ set (fig.width = 7 , # 7" width fig.asp = 0.618 , # the golden ratio fig.retina = 3 , # dpi multiplier for displaying HTML output on retina fig.align = "center" , # center align figures dpi = 300 # higher dpi, sharper image = TRUE ``` ## Plan out layer mappings - geom_point:``` - Creating a scatter plot of the seismic events (depth and duration-magnitude) over time ``` - geom_line: - Creating a scatter plot of the seismic events (depth and duration-magnitude) over time- geom_bar: - View the frequency of the seismic events aggregated per year or year, month```{r} <- read.csv (here ("data" ,"vesuvius.csv" ))<- vesuvius_data |> mutate (date = str_replace_all (time, "T" , " " ),date = str_replace_all (date, "Z" , "" ),date = as.POSIXct (date, format = "%Y-%m-%d %H:%M:%S" )<- vesuvius_data %>% filter (! is.na (duration_magnitude_md) & ! is.na (depth_km))glimpse (vesuvius_data)glimpse (vesuvius_data_filtered)``` ```{r} $ year <- format (vesuvius_data_filtered$ date, "%Y" )$ month <- format (vesuvius_data_filtered$ date, "%m" )$ year <- as.numeric (as.character (vesuvius_data_filtered$ year))$ year <- as.numeric (as.character (vesuvius_data_filtered$ year))``` ```{r} <- vesuvius_data_filtered %>% group_by (year) %>% summarise (mean_duration_magnitude_md = mean (duration_magnitude_md), mean_depth_km = mean (depth_km))``` ```{r} # aggregated_data$date <- paste0(aggregated_data$year, " ", aggregated_data$month) glimpse (aggregated_data)``` ```{r} <- ggplot (aggregated_data, aes (x = year, y = mean_duration_magnitude_md)) + geom_jitter () + geom_smooth (method = "loess" , se = FALSE ) + scale_x_continuous (breaks = pretty (aggregated_data$ year, n = 10 ))+ labs (title = "Duration Magnitude averaged per year" , caption = "Source: Seismic data for Mount Vesuvius \n TidyTuesday 2025-05-13 \n Graphic: Cedric Destin" ,x = "Year" , y = "Duration Magnitude" )if (save_figs) {ggsave ("images/duration_magnitude_by_year.png" , plot_1)``` ```{r} <- ggplot (aggregated_data, aes (x = year, y = mean_depth_km)) + geom_jitter () + geom_smooth (method = "loess" , se = FALSE ) + scale_x_continuous (breaks = pretty (aggregated_data$ year, n = 10 )) + labs (title = "Depth averaged per year" , caption = "Source: Seismic data for Mount Vesuvius \n TidyTuesday 2025-05-13 \n Graphic: Cedric Destin" ,x = "Year" , y = "Depth Km" )if (save_figs) {ggsave ("images/depth_by_year.png" , plot_2)``` ```{r} <- vesuvius_data_filtered %>% group_by (year, month) %>% summarise (mean_duration_magnitude_md = mean (duration_magnitude_md), mean_depth_km = mean (depth_km))``` ```{r} $ date <- paste0 (aggregated_data_2$ year, " " , aggregated_data_2$ month)glimpse (aggregated_data)``` ```{r} $ date <- as.Date (with (aggregated_data_2, paste (year, month, "01" , sep = "-" )))<- ggplot (aggregated_data_2, aes (x = date, y = mean_depth_km)) + geom_jitter () + geom_smooth (method = "loess" , se = TRUE ) + scale_x_date (date_labels = "%Y" , date_breaks = "12 month" ) + labs (title = "Depth averaged per year & month" ,y = "Depth Km" , x = NULL ) + theme (axis.text.x = element_blank (), axis.ticks.x = element_blank ())if (save_figs) {ggsave ("images/depth_by_year_month.png" , plot_3)``` ```{r} <- ggplot (aggregated_data_2, aes (x = date, y = mean_duration_magnitude_md)) + geom_jitter () + geom_smooth (method = "loess" , se = TRUE ) + scale_x_date (date_labels = "%Y" , date_breaks = "12 month" ) + labs (title = "Duration Magnitude averaged per year & month" , caption = "Source: Seismic data for Mount Vesuvius \n TidyTuesday 2025-05-13 \n Graphic: Cedric Destin" ,x = "Year" , y = "Duration Magnitude" )if (save_figs) {ggsave ("images/duration_magnitude_by_year_month.png" , plot_4)``` ```{r} <- grid.arrange (plot_3, plot_4, nrow = 2 )if (save_figs) {ggsave ("images/by_year_month.png" , plot_3_4)``` ```{r} <- ggplot (filter (vesuvius_data_filtered, year >= 2013 ), aes (x = depth_km)) + geom_histogram (fill = "skyblue" , color = "black" , bins = 100 ) + labs (title = "Histogram of Seismic Depth" , x = "Depth (km)" , y = "Count of Seismic Events" , caption = "Source: Seismic data for Mount Vesuvius \n TidyTuesday 2025-05-13 \n Graphic: Cedric Destin" )+ facet_wrap (~ year)if (save_figs) {ggsave ("images/depth_hist.png" , plot_5)<- ggplot (filter (vesuvius_data_filtered, year >= 2013 ), aes (x = duration_magnitude_md)) + geom_histogram (fill = "skyblue" , color = "black" , bins = 100 ) + labs (title = "Histogram of Seismic Magnitude Duration" , x = "Duration Magnitude" , y = "Count of Seismic Events" , caption = "Source: Seismic data for Mount Vesuvius \n TidyTuesday 2025-05-13 \n Graphic: Cedric Destin" ) + facet_wrap (~ year)if (save_figs) {ggsave ("images/magnitude_duration_hist.png" , plot_6)``` ```{r include=FALSE} # p1 <- ggplot(filter(vesuvius_data_filtered, year == 2013), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p2 <- ggplot(filter(vesuvius_data_filtered, year == 2014), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p3 <- ggplot(filter(vesuvius_data_filtered, year == 2015), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p4 <- ggplot(filter(vesuvius_data_filtered, year == 2016), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p5 <- ggplot(filter(vesuvius_data_filtered, year == 2017), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p6 <- ggplot(filter(vesuvius_data_filtered, year == 2018), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p7 <- ggplot(filter(vesuvius_data_filtered, year == 2019), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p8 <- ggplot(filter(vesuvius_data_filtered, year == 2020), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p9 <- ggplot(filter(vesuvius_data_filtered, year == 2021), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p10 <- ggplot(filter(vesuvius_data_filtered, year == 2022), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p11 <- ggplot(filter(vesuvius_data_filtered, year == 2023), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p12 <- ggplot(filter(vesuvius_data_filtered, year == 2024), aes(x = depth_km)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, ncol = 3) ``` ```{r include=FALSE} # p1 <- ggplot(filter(vesuvius_data_filtered, year == 2013), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p2 <- ggplot(filter(vesuvius_data_filtered, year == 2014), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p3 <- ggplot(filter(vesuvius_data_filtered, year == 2015), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p4 <- ggplot(filter(vesuvius_data_filtered, year == 2016), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p5 <- ggplot(filter(vesuvius_data_filtered, year == 2017), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p6 <- ggplot(filter(vesuvius_data_filtered, year == 2018), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p7 <- ggplot(filter(vesuvius_data_filtered, year == 2019), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p8 <- ggplot(filter(vesuvius_data_filtered, year == 2020), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p9 <- ggplot(filter(vesuvius_data_filtered, year == 2021), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p10 <- ggplot(filter(vesuvius_data_filtered, year == 2022), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p11 <- ggplot(filter(vesuvius_data_filtered, year == 2023), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # p12 <- ggplot(filter(vesuvius_data_filtered, year == 2024), aes(x = duration_magnitude_md)) + # geom_histogram(fill = "skyblue", color = "black", bins = 100) + # labs(title = "Histogram of your_column", x = "your_column", y = "Count") # grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, ncol = 3) ``` ```{r include=FALSE} # ggplot(vesuvius_data_filtered, aes(x = date, y = duration_magnitude_md)) + # geom_point(position = position_dodge(width = 0.5)) + # stat_dist_halfeye( # aes(dist = dist_normal(mu = AME, sigma = SE)), # color = "#0072B2", # fill = "#0072B288", # alpha = NA, # point_size = 4 # ) + # facet_grid(rows = vars(date), # # labeller = labeller(response = label_wrap_gen(15), # # explanatory = label_wrap_gen(15)), # space = "free_y", # scales = "free_y") # labs(x = "Mean likert score \n (Error bars range from 10th to 90th percentile)", y = "") + # theme_minimal() + # theme( # strip.background = element_rect(fill = "gray90", color = "lightgray"), # strip.text.x = element_text(angle = 0), # strip.text.y = element_text(angle = 0) # ) ``` ```{r include=FALSE} # Replace 'aggregated_data', 'x_column', and 'y_column' with your data and columns <- ggpairs (aggregated_data_2, columns = c ("mean_depth_km" , "mean_duration_magnitude_md" ),upper = list (continuous = wrap ("points" , alpha = 0.5 )),lower = list (continuous = wrap ("cor" , size = 3 )),title = "Pairwise Scatterplot Matrix of Depth and Duration Magnitude" )if (save_figs) {ggsave ("images/pairwise_plot.png" , plot_7)``` ```{r} <- vesuvius_data_filtered %>% group_by (year) %>% summarize (avg_depth_km = mean (depth_km, na.rm = TRUE ),sd_depth_km = sd (depth_km, na.rm = TRUE ),max_depth_km = max (depth_km, na.rm = TRUE ),min_depth_km = min (depth_km, na.rm = TRUE ),avg_duration_magnitude_md = mean (duration_magnitude_md, na.rm = TRUE ),sd_duration_magnitude_md = sd (duration_magnitude_md, na.rm = TRUE ),max_duration_magnitude_md = max (duration_magnitude_md, na.rm = TRUE ),min_duration_magnitude_md = min (duration_magnitude_md, na.rm = TRUE )``` ```{r} <- ggplot (vesuvius_data_summary) + aes (y = avg_depth_km, x = year) + # geom_vline(xintercept = 0, color = "gray20") + stat_dist_halfeye (aes (dist = dist_normal (mu = avg_depth_km, sigma = sd_depth_km)),color = "#0072B2" ,fill = "#0072B288" ,alpha = NA ,point_size = 4 if (save_figs) {ggsave ("images/duration_magnitude_hist.png" , plot_8)``` ```{r} <- ggplot (vesuvius_data_summary) + aes (y = avg_duration_magnitude_md, x = year) + geom_point () + geom_linerange (aes (ymin = min_duration_magnitude_md, ymax = max_duration_magnitude_md)) + geom_errorbar (aes (ymin = min_duration_magnitude_md, ymax = max_duration_magnitude_md), width = 0.2 ) + geom_smooth (method = "loess" , se = FALSE ) + labs (title = "Average Duration Magnitude by Year" ,caption = "Source: Seismic data for Mount Vesuvius \n TidyTuesday 2025-05-13 \n Graphic: Cedric Destin" ,x = "Year" , y = "Average Duration Magnitude" )if (save_figs) {ggsave ("images/duration_magnitude_error.png" , plot_9)``` ```{r} <- ggplot (vesuvius_data_summary) + aes (y = avg_depth_km, x = year) + geom_point () + geom_linerange (aes (ymin = min_depth_km, ymax = max_depth_km)) + geom_errorbar (aes (ymin = min_depth_km, ymax = max_depth_km), width = 0.2 ) + geom_smooth (method = "loess" , se = FALSE ) + labs (title = "Average Depth by Year" ,caption = "Source: Seismic data for Mount Vesuvius \n TidyTuesday 2025-05-13 \n Graphic: Cedric Destin" ,x = "Year" , y = "Average Depth (km)" )if (save_figs) {ggsave ("images/depth_error.png" , plot_10)``` ```{r} head (vesuvius_data_filtered)$ binned_duration <- cut ($ duration_magnitude_md,breaks = seq (floor (min (vesuvius_data_filtered$ duration_magnitude_md, na.rm = TRUE )),ceiling (max (vesuvius_data_filtered$ duration_magnitude_md, na.rm = TRUE )),by = 0.25 include.lowest = TRUE $ bin_midpoint <- as.numeric (sub (" \\ ((.+),(.+) \\ ]" , " \\ 1" , levels (vesuvius_data_filtered$ binned_duration)))[vesuvius_data_filtered$ binned_duration] + 0.25 ``` ```{r} <- vesuvius_data_filtered %>% group_by (year, binned_duration, bin_midpoint) %>% summarise (mean_duration_magnitude_md = mean (duration_magnitude_md), count_duration_magnitude_md = n ())``` ```{r} glimpse (aggregated_data_3)``` ```{r} <- ggplot () + # Layer 1: vesuvius_data_summary (blue) geom_point (data = vesuvius_data_summary, aes (x = year, y = avg_depth_km, color = "Average Depth" ), size = 3 , shape = 1 ) + geom_linerange (data = vesuvius_data_summary, aes (x = year, ymin = min_depth_km, ymax = max_depth_km, color = "Depth Range" )) + geom_errorbar (data = vesuvius_data_summary, aes (x = year, ymin = min_depth_km, ymax = max_depth_km, color = "Depth Range" ), width = 0.2 ) + geom_smooth (data = vesuvius_data_summary, aes (x = year, y = avg_depth_km, color = "Trend" ), method = "loess" , se = FALSE ) + # Layer 2: aggregated_data_3 (red, size mapped) geom_point (data = aggregated_data_3, aes (x = year, y = bin_midpoint, size = count_duration_magnitude_md, color = "Depth Bin" ), alpha = 0.6 ) + scale_color_manual (name = "Layer" ,values = c ("Average Depth" = "blue" ,"Depth Range" = "darkblue" ,"Trend" = "skyblue" ,"Depth Bin" = "red" + scale_size_continuous (name = "Count Seismic Events" ,guide = guide_legend (override.aes = list (color = "red" ))+ labs (size = "Count Seismic Events" , y = "Depth Km" , x = "Year" , caption = "Source: Seismic data for Mount Vesuvius \n TidyTuesday 2025-05-13 \n Graphic: Cedric Destin" ) + scale_x_continuous (breaks = pretty (aggregated_data$ year, n = 10 ))if (save_figs) {ggsave ("images/depth_error_binned.png" , plot_11)``` ```{r} # Example overlay plot with legends for each layer <- ggplot () + # Layer 1: vesuvius_data_summary (blue) geom_point (data = vesuvius_data_summary, aes (x = year, y = avg_duration_magnitude_md, color = "Average Duration Magnitude" , size = 50 ), shape = 1 ) + geom_linerange (data = vesuvius_data_summary, aes (x = year, ymin = min_duration_magnitude_md, ymax = max_duration_magnitude_md, color = "Duration Magnitude Range" )) + geom_errorbar (data = vesuvius_data_summary, aes (x = year, ymin = min_duration_magnitude_md, ymax = max_duration_magnitude_md, color = "Duration Magnitude Range" ), width = 0.2 ) + geom_smooth (data = vesuvius_data_summary, aes (x = year, y = avg_duration_magnitude_md, color = "Trend" ), method = "loess" , se = FALSE ) + # Layer 2: aggregated_data_3 (red, size mapped) geom_point (data = aggregated_data_3, aes (x = year, y = bin_midpoint, size = count_duration_magnitude_md, color = "Duration Bin" ), alpha = 0.6 ) + scale_color_manual (name = "Layer" ,values = c ("Average Duration Magnitude" = "blue" ,"Duration Magnitude Range" = "darkblue" ,"Trend" = "skyblue" ,"Duration Bin" = "red" + scale_size_continuous (name = "Count Duration Magnitude" ,guide = guide_legend (override.aes = list (color = "red" ))+ labs (title = "Duration Magnitude of Seismic Events at Mount Vesuvius Over Time" ,size = "Count Duration Magnitude" , y = "Duration Magnitude" , x = "Year" ,caption = "Source: Seismic data for Mount Vesuvius \n TidyTuesday 2025-05-13 \n Graphic: Cedric Destin" )+ scale_x_continuous (breaks = pretty (aggregated_data$ year, n = 10 ))if (save_figs) {ggsave ("images/duration_magnitude_error_binned.png" , plot_12)``` ```{r} <- (filter (vesuvius_data, ! is.na (longitude) & ! is.na (latitude)))<- vesuvius_data_filter_2 %>% filter (! is.na (longitude)) %>% select (x = longitude, y = latitude) %>% as.data.frame () %>% get_elev_raster (locations = ., z = 13 ,prj = 4326 ``` ```{r} <- terra:: as.data.frame (vesuvius_hm, xy = T) %>% na.omit () %>% rename (elevation = 3 )``` ```{r} # Create a diagonalsection and get the elevation <- 14.34 <- 40.785 <- 14.50 <- 40.85 # Tolerance for getting the points on the line <- 5e-6 # Calculate the cross product to find points on the line <- vesuvius_hm_df %>% mutate (cross_product = (y - p1_y) * (p2_x - p1_x) - (x - p1_x) * (p2_y - p1_y)%>% filter (abs (cross_product) < epsilon & >= (min (p1_x, p2_x) - epsilon) & x <= (max (p1_x, p2_x) + epsilon) & >= (min (p1_y, p2_y) - epsilon) & y <= (max (p1_y, p2_y) + epsilon)``` ```{r} $ time <- as.POSIXct (vesuvius_data_filter_2$ time)<- min (as.numeric (vesuvius_data_summary$ year), na.rm = TRUE )<- max (as.numeric (vesuvius_data_filter_2$ year), na.rm = TRUE )``` ```{r} <- v_section %>% mutate (x_rescaled_to_date = rescale (x, to = c (date_min, date_max))``` ```{r} <- v_section %>% mutate (x_rescaled_to_date = as.POSIXct (x_rescaled_to_date, origin = "1970-01-01" )``` ```{r} <- ggplot () + geom_segment (data = v_section, aes (x = x_rescaled_to_date, xend = x_rescaled_to_date, y = min (elevation / 1000 ), yend = elevation / 1000 ), color = "#eed9c4" ) + # Layer 1: vesuvius_data_summary (blue) geom_point (data = vesuvius_data_summary, aes (x = year, y = avg_depth_km*- 1 , color = "Average Depth" , size = 50 ), shape = 1 ) + geom_linerange (data = vesuvius_data_summary, aes (x = year, ymin = min_depth_km*- 1 , ymax = max_depth_km*- 1 , color = "Depth Range" )) + geom_errorbar (data = vesuvius_data_summary, aes (x = year, ymin = min_depth_km*- 1 , ymax = max_depth_km*- 1 , color = "Depth Range" ), width = 0.2 ) + geom_smooth (data = vesuvius_data_summary, aes (x = year, y = avg_depth_km*- 1 , color = "Trend" ), method = "loess" , se = FALSE ) + # Layer 2: aggregated_data_3 (red, size mapped) geom_point (data = aggregated_data_3, aes (x = year, y = bin_midpoint*- 1 , size = count_duration_magnitude_md, color = "Depth Bin" ), alpha = 0.6 ) + scale_color_manual (name = "Layer" ,values = c ("Average Depth" = "blue" ,"Depth Range" = "darkblue" ,"Trend" = "skyblue" ,"Depth Bin" = "red" + scale_size_continuous (name = "Count Seismic events" ,guide = guide_legend (override.aes = list (color = "red" ))+ labs (title = "Depth of Seismic Events at Mount Vesuvius Over Time" ,size = "Count Seismic Events" , y = "Depth Km" , x = "Year" , caption = "Source: Seismic data for Mount Vesuvius \n TidyTuesday 2025-05-13 \n Graphic: Cedric Destin" )+ scale_x_continuous (breaks = pretty (aggregated_data$ year, n = 10 ))if (save_figs) {ggsave ("images/depth_error_binned_mount.png" , plot_13)``` ## Plan out data wrangling methods - time: - Replace "T" with " " - Remove "Z" - Convert into a date variable - Extract the years - Extract the months
